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Abstract 

This paper presents an approach tor the generation ot closed mani- 
fold surface triangulations trom CAD geometry. CAD parts and 
assemblies are used in their native format, without translation, and a 
part’s native geometry engine is accessed through a modeler-inde- 
pendent application programming interlace (API). In seeking a 
robust and fully automated procedure, the algorithm is based on a 
new physical space manitold triangulation technique which was 
developed to avoid robustness issues associated with poorly condi- 
tioned mappings. In addition, this approach avoids the usual ambi- 
guities associated with floating-point predicate evaluation on 
constructed coordinate geometry in a mapped space. The technique 
is incremental, so that each new site improves the triangulation by 
some well defined quality measure. Sites are inserted using a variety 
of priority queues to ensure that new insertions will address the 
worst triangles first. As a result ot this strategy, the algorithm will 
return its “best” mesh for a given (prespecified) number of sites. 
Alternatively, the algorithm may be allowed to terminate naturally 
after achieving a prespecified measure ot mesh quality. The result- 
ing triangulations are “CFD-ready” in that: (l) Edges match the 
underlying part model to within a specified tolerance. (2) Triangles 
on disjoint surfaces in close proximity have matching length-scales. 
(3) The algorithm produces a triangulation such that no angle is less 
than a given angle bound, cx, or greater than k - 2oc. This result also 
sets bounds on the maximum vertex degree, triangle aspect-ratio 
and maximum stretching rate tor the triangulation. In addition to the 
output triangulations for a variety of CAD parts, the discussion pre- 
sents related theoretical results which assert the existence ot such an 
angle bound, and demonstrate that maximum bounds of between 
25° and 30° may be achieved in practice. 

1. Introduction 

Mesh generation has long been recognized as a bottleneck in 
the CFD process. [ 11 The last decade has witnessed a myriad 
of international and domestic conferences and symposiums 
aimed at focusing research on this impediment. Unstructured, 
hybrid, and Cartesian mesh methods are all aimed at simpli- 
fying the mesh generation task for complex configurations. 
The success of these approaches is well represented in the lit- 
erature and with an appropriate initial surface triangulation, 
the volume mesh generation can generally be accomplished 
in a relatively automated fashion in minutes-to-hours on an 
engineering workstation.^ As taster processors, with bet- 
ter access to memory, continue to shrink the wall-clock time 
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required for both mesh generation and flow solution, the 
man-hour intensive task of generating an initial surface trian- 
gulation from a CAD geometry promises to become an ever 
larger fraction of the bottleneck. Additionally, it the user 
must be involved in the extraction of surface data from CAD, 
then mesh adaptation - which involves enriching the discreti- 
zation on the body surface - will remain an elusive goal. 

Historically, surface discretization has been one of the least 
automated steps in the numerical simulation cycle, and tor 
good reason. Due to its dependence on implicitly defined sur- 
faces and curves, CAD data is by its nature imprecise. Vari- 
ous geometry engines typically demonstrate discrepancies in 
their interpretations of the same entities. As a result, “repair 
of CAD surfaces has become an area of substantial 
researchJ 8 ^ This problem is exacerbated when CAD mod- 
els are output in many of the standard formats, since such 
files frequently do not include important topological and con- 
struction information along with the entity geometry. 

In addition to the vagaries of the data and modeling systems, 
many surface triangulation schemes treat surfaces individu- 
ally. Thus, the length scale of the triangulation on a particular 
component may not be commensurate with that of other com- 
ponents which are in close proximity on the full configura- 
tion. In many interactive triangulation systems, it becomes 
the user’s role to identity such situations and enforce length 
scale compatibility between nearby components. 

In response to these and other requirements for user assis- 
tance, many in the research and industrial CFD communities 
have adopted an interactive paradigm for surface mesh gener- 
ation. The commercial unstructured mesh generators in Refs. 
[7], [10] and [11] all interact with CAD data through files 
which have been translated from their CAD native environ- 
ment to some standardized format (namely IGES , 
STEP 1131 or STlJ 14] ). 1 

This paper adopts an alternative paradigm. The approach 
interfaces with the CAD system in a dynamic manner via 
calls to CAD native routines. By accessing the model in its 
native environment, this approach avoids translation to a for- 
mat which can deplete the model ot topological information. 
This is important since it avoids the consistency conflicts that 


1 Recent releases of some ot this software now supports direct 
interfaces which do read parts in their native formats, however, this 
practice is not the norm. 
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can occur when two different geometry engines attempt to 
infer topological information from imprecise data. 

To avoid placing CAD specific calls in the software, we 
access CAD function calls through an Application Program- 
ming Interface (API) known as CAPRI. 1 2 ' 51 This library pro- 
sents a standardized interface to the application program ior 
various CAD systems. 1 CAPRI supports a variety of opera- 
tions like truth testing, geometry construction, and entity 
queries. 

Maintaining the consistency of the models by direct manipu- 
lation of CAD parts and assemblies is the first step that this 
work takes toward building a robust method for surface tri an- 
gulation. A basic premise of this approach is that we resolve 
consistency problems on as simple a model as possible, and 
maintain this consistency as the triangulation evolves and 
becomes more complex. 

1.1. Abstract Geometric Structures 

Following the approach of YapJ 16 ^ our approach toward gen- 
erating a robust geometric algorithm contends that a geomet- 
ric structure, D, consists of four elements: - 

D = (c, X, O(z), z ) (1) 

Where the graph, G = (V, E) is a directed set of vertices, V, 
and edges, E. X is a function describing the index labels of 
the graph. O is a geometric operator which represents the 
consistency predicates for the connectivity and is a function 
of the actual coordinates z. G represents a tessellation ol the 
vertices and is therefore purely combinatorial. A structure is 
said to be consistent if the predicate O(z) holds. 

As an example, a 2-D Delaunay triangulation algorithm usu- 
ally makes use of an inCircle predicate^ 17 ], ^>i n cirde which 
establishes G by insisting that the circumcircle of the triangle 
A a i, c can contain no other vertex in the graph. If such a pred- 
icate holds for every triangle in G, then this instance of the 
geometric structure D is said to be consistent. 

1.2. Robustness 

Consistent CAD Models 

This interpretation offers direct insights into the formulation 
of a robust algorithm for creating triangulations of CAD vol- 
umes. The rational B-splines used to describe surfaces in 
most CAD systems are implicitly defined for physical space 
coordinates of the geometry. Therefore, the constructors for 
vertex geometry generally require a Newton solve carried to 
some internal tolerance. Since the results of this construction 
will be subject to both tolerance and round-off error, the sys- 
tem may then “nudge” the constructed point, z f , to some 
nearby exactly representable location (on an integer grid, for 
example). If the geometry engine’s predicate for determining 
if Zj is on some surface, S is O ( z ;> ^), is consistent then it 

1 onSurf*- 


1 CAPRI currently supports ProEngineer™ ParaSolids™ and 
IDEAS™ with CATI A™ support in beta test. 

2 Ref.[l6] actually writes eq.(l) as D- (G, X, O(z), l) where / isjt 

mapping from the input parameters c to z, /; z->c = (< , • 


will return “trup” when later queried if z, lies on the surface. 
However, if & onSurf now represents some user-defined 
predicate which may be ignorant of the systems construction 
rules, then it is very unlikely to return consistent results. 

The CAPRI API allows us to access a subset of the CAD 
geometry engine’s constructors, queries and predicates in 
order for the algorithm in our application to maintain a con- 
sistent representation of the model. Implementationally, we 
adopt a multi-threaded programming approach which starts 
the geometry engine on its own thread in order to respond to 
queries from the main triangulation thread. 

Physical Space Triangulation 

A variety of existing surface meshing techniques adopt a 
mapped- space approach for generating surface triangula- 
tions. In this approach, the surface and its bounding curves 
are triangulated in a 2-D parameter space, which may or may 
not have some additional scaling imposed. Ref. f 18] provides 
a mathematical description of the Non-Uniform Rational B- 
Spline (NURBS) surfaces typically used in CAD systems. 
Here we note only that an iterative method is required to 
solve for the physical space coordinates of a position speci- 
fied on the surface in the parameter space. This process 
involves division of two (generally) high-order polynomials, 
and it therefore incurs both floating-point round-off error and 
error associated with tolerancing for the convergence of the 
iterative solve. As a result, computed coordinates in the 
mapped space are necessarily noisy and cannot be considered 
exact values. Two consequences of this approach are: 

1. Since the error bounds on the input are unknown, evalu- 
ation of the triangulation predicate (e.g. $>i n circle ) are 
unlikely to robustly produce consistent results (see Fig. 

19 Ref.[ 17], also [19] and [20]). 

2. The polynomial basis for the NURB may be high-order, 
and therefore small errors in parameter space may pro- 
duce dramatic results in physical space - even within the 
subspace for which the surface is defined. The likeli- 
hood of encountering poorly conditioned mappings is 
the primary reason that CAD repair software generally 
attempts to re-normalize the NURBS surface and recast 
it using basis polynomials with as low an order as possi- 
ble^. 

These two observations have led the current work to focus on 
physical space triangulation techniques. In this approach, we 
construct a manifold triangulation on the surface, and evalu- 
ate the triangulation predicates in New sites are con- 
structed by the CAD geometry engine and, since this output 
is consistent with the system’s internal predicates, it is con- 
sidered exact by the external predicates of the triangulation 
algorithm. The presentation in §3 places emphasis on mini- 
mizing the floating-point error in predicate evaluation and 
notes the computational requirements needed to make this 
evaluation exact. 

2. The CAPRI API 

Our basic approach is to take a crude manifold triangulation 
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of each closed volume in the CAD assembly, and improve it 
until it satisfies a preset measure of mesh quality, or produces 
a preset number of triangles. A variety of mesh quality mea- 
sures may be defined within this framework, and this prelim- 
inary investigation focuses on three such criteria. (1) The 
mesh must be free from small angles (sliver triangles). ( 2 ) 
Edges in the triangulation must not deviate from the underly- 
ing model by more than a prescribed tolerance. (3) Length 
scales on nearby (possibly disjoint) surfaces must be com- 
mensurate. 


versed in a counterclockwise circuit when viewed from a 
point outside the solid. This convention permits holes in a 
surface to be described by a clockwise loop. In the figure, 
cad_face,//, consists of loops lj and l 2 , each of which is com- 
posed of cad_edge entities. The edge ordering of l 2 indicates 
that it is clockwise, and therefore describes a hole in /y. 
Edges and faces have an underlying parameterization, and 
while points may be queried for their parametric values 
(u, v), details of this ruling are not exposed to the application 
program. 


2.1. CAPRI Volumes 


2.2. Initial Manifold Triangulation 


As discussed in the introduction, CAD entities are accessed 
through the CAPRI programming interfaced 131 CAPRI pro- 
vides a layer of indirection such that CAD system specific 
data may be accessed by an application program using CAD 
system neutral function calls. Figure 1 presents an abstract 
view of the entities that CAPRI provides. A cadjnode is the 
lowest dimensional entity and corresponds to a point in 3- 
space. A cad_edge has a cad_node at both ends. Each is 
directed from its origin, 0, to its destination, D. Cad_edges 
are not assumed to be simplicial and may follow a general 
curve in space (see e$ and e 5 in Fig.l). Each edge is con- 
nected to two cad_face entities. In general, these faces are 
composed of several loops and are not assumed planar, since 
they follow the underlying parameterization of the surface. A 
cad_face is composed of one or many loops , which are col- 
lections of oriented edges. Loops are oriented such that the 
surface of the cad_face lies inside them when they are tra- 

2 .a 

l 


A central theme in the present approach is the maintenance of 
a closed volume throughout the procedure. CAPRI returns a 
simplicial decomposition of each of the m cad_face entities, 
S h where 5- € {5,, S 2 , ..., S m } . Each of these triangulations 
are manifold within their respective cad_cdges. In addition, 
an indexing function, X c , is returned for each cad_volume. 
Therefore a simplicial, manifold representation of each 
cad_volume, S c may be constructed by taking the union of 
the decompositions of all the cad_face entities of a volume, 
subject to the indexing X c . 

S c = ^i s i V/ e { 1, ..., m} (2) 

Figure 2 displays an example of this initial triangulation for a 
simple part. The manufacturing die shown has 14 cad_face 
entities and the initial triangulation, has 270 triangles. As 
is typical, this triangulation is quite irregular, and planar 
regions are decomposed into as few triangles as possible. 




Figure 2. Initial dosed, manifold, surface triangulation, S ( . t of a CAD model for a manufacturing die. Underlying CAD 
entities exposed to application program are labeled in the frame on the right. 
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Extremely high aspect ratio triangles are also common in 
such boundary triangulations. Figure 2.b labels selected CAD 
entities on this triangulation. Notice that although some 
cad_face sites may be present, this initial triangulation is 
essentially a boundary triangulation and the number of trian- 
gles is proportional to the number of cad„edges. 

Despite the poor quality of the triangulation, the structure in 
Fig. 2 has several desirable properties. Namely, it is manifold, 
oriented and closed. We wish to improve this triangulation by 
adding sites on both the cad_edge and cad_face entities and 
by enforcing an external predicate governing the type of tri- 
angulation. 

3. Mesh Improvement 


Our approach for manifold surface triangulation traces its 
roots to work on quality triangulations of Planar Straight 
Line Graphs (PSLGs) [21][22][23] and related work on quality 
triangulations of manilold surfaces J 24 ^ Work in this field 
began with the efforts of Ref.[21] which presented an algo- 
rithm with both shape and size guarantees. The resulting 
meshes were size-optimal and had no triangle with an aspect 
ratio greater than 5. In this context, the aspect ratio AR , of a 
triangle, is defined as the length of the longest edge divided 
by that of the shortest one. One can show that if a is the 
smallest angle of a triangle, then 


1 


< AR < i 


(3) 


| sin oc| |sina[ 

Therefore a is frequently used to describe the quality of a 
given triangle. 


Before presenting the manifold triangulation technique, this 
section first recounts a related algorithm for quality triangula- 
tion for PSLGs from Ref.[23]. It then present a fundamen- 
tally similar algorithm for triangulating curved surfaces and 
note which aspects of the PLSG method have been relaxed in 
the extension. 


3.1. Quality Triangulation of PSLGs 

While the manifold surface triangulation technique of Rup- 
pert [24] and the PSLG method of Chew [23] are similar in 
many respects, our manifold technique follows Ruppert s 
approach more closely. Section 3.2 addresses some of the 
reasoning behind this choice. 

An essential feature of the algorithm is the notion of an 
encroached constraining edge. As illustrated in Figure 3, a 
constraining edge, OD , is said to be encroached upon if any 
other site (visible to OD ) lies within the diametral circle of 
the edge. 

If one recalls that the circumcenter of a right triangle falls on 
the hypotenuse, then its easy to show that for a triangulation 
which is Delaunay or locally maxmin, a predicate for 
encroachment may be formulated as a vector dot product. 
Thus for similarly sized 1 floating-point data with p-bit 
significands, this predicate can be computed exactly in a reg- 
ister with a 2p-bit significant 201 . In a practical sense, this 
implies that as long as the edges are small with respect to 
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Figure 3. Constraining edge OD and its diametral circle. The 
edge is “ encr oached” if any site, p , falls within the diametral 
circle of OD 


their distance from the origin, this predicate will be exact if 
computed in double-precision, using single-precision data. 

In the algorithms, Encroached^’) denotes application of this 
predicate to an edge, e. The (~) superscript is used to remind 
us that since this predicate is part of our triangulation algo- 
rithm, it is not native to the CAD system. 

The presentation of Ref.[23] recovers the constraining edges 
of the triangulation as the algorithm advances. To clarify its 
relation with our manifold triangulation technique, we recast 
the original algorithm assuming that it begins with a con- 
strained boundary triangulation of the input vertices, V, of the 
constraint edges. Furthermore, this initial triangulation is 
assumed to be the constrained Delaunay triangulation of the 
input sites, CL/JfV). 

The algorithm is quite elegant in that it consists of only two 
major operations: 

1 . Split a constrained edge: Add a site to the mid-point of a 
constraining edge, and replace the original edge with the 
two new edges in the constraint list. 

2. Split a triangle: Add a site to the circumcenter of a trian- 
gle, /. Note that if the triangle is obtuse, this site will not 
fall within t. 

Algorithm Q: Quality triangulation of a PLSG 

Input: Planar Straight Line Graph, X, with input vertices, V m . 

Target angle, a. 

Output: CmV nll( ) with all angles > a. 

Initialize: Compute ODOR Build minimum angle priority 
queue, PQ mm with t P g denoting triangle at head of 
queue, having min angle 0 P q 

1. Apply <S> Encroached to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CJXRV), Update PQ mm - 

2. While < OC) { 

2. a Let p be the circumcenter of t P Q. 

2.b If (p encroaches any constraining edge, e) 

2.c Split constrained edge. Update C‘M(V), Update PQ m i tv 

2. d Else Split triangle t P Q. 

Add p to L Update CM XV), Update PQ min - 

1 

3. Output CM[V). 


1. The qualifier “similarly sized” is necessary to guard against the 
case where a coordinates of one point is less than half that of 
another. Extended precision would be required in such a case. 
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While simple, Ref.[23] proves that Alg.Q produces triangula- 
tions with the following desirable properties. 

1 . Quality: An angle bound 20.7° is guaranteed, and values as 
high as 30° may be achieved in practice. 

2. Output Size: The size of the output triangulation is within 
a constant of the optimal number of triangles required to sat- 
isfy the angle criteria. 

3. Size Optimality: Small input constraint edges are sur- 
rounded by proportionally small triangles. Nearby triangles 
have similar sizes, and the size variation of triangles in the 
mesh is proportional to the distance between them. 

3.2. Difficulties on Curved Surfaces 

In Ref. [24], Chew presents a quality triangulation technique 
which is closely related to Alg.Q. This work raises a number 
of difficulties associated with the extension of the PSLG 
method to curved surfaces. 

Alg.Q has two salient aspects. (1) The triangulation is con- 
strained Delaunay. (2) New sites are added at circumcenters. 
In making these observations, Chew notes that a straightfor- 
ward definition of a circle on a surface is the loci of points on 
the surface which are equidistant from another point on the 
surface, where all distances are measured using the geodesic 
distance along the surface. While straightforward, this defini- 
tion is problematic. Distances along the surface must be mea- 
sured in physical space, and will therefore be expensive to 
compute on NURBS surfaces. In addition, due to the inherent 
error in finding the coordinates of a point on such a surface, 
robust predicates based upon this definition will be difficult 
to formulate. Finally, Chew notes that this definition has less 
subtle and non-intuitive consequences. A circle whose center 
lies near the base of a sharp spire, for example, may reach 
completely around the spire without also including the tip of 
the spire. 

To circumvent such difficulties, the method in [24] makes use 
of an alternative definition of a circle. In the plane, the three 
vertices of a triangle define a unique circle. In 3-dimensions, 
however, an infinite family of spheres may be passed through 
those three points. Connecting the line through the centers of 
this family of spheres and intersecting this line with the sur- 
face identifies a particular sphere in this family. The circum- 
center of the triangle may be defined to be the loci of points 
at the intersection of this particular sphere and the surface. 
Once this sphere is found, then the inCircle triangulation 
predicate may be evaluated by simply computing distances in 
three dimensions as a vector magnitude. 

Despite the effort, problems still exist with this predicate. ( 1 ) 
If the triangulation does not resolve the underlying surface 
closely enough, the line of circumsphere centers will not nec- 
essarily intersect the surface. Alternatively, it may also inter- 
sect the surface in multiple locations places. (2) Computing 
the intersection of this line with the surface will require an 
iterative method and will therefore be computationally inten- 


sive. (3) Once this intersection is successfully located on the 
surface, one must determine which triangle the point falls 
inside. Since the triangulation only matches the surface at the 
vertices, ambiguous situations may arise when two triangles 
claim ownership of the same site (near a ridge, for example). 
(4) On a curved surface, the circumcenters of two triangles 
with a shared edge may not be consistent. Specifically, when 
(J > inCircle tests one °f the triangles against the opposite ver- 
tex of the other, the results may not agree when the roles of 
the triangles are reversed. Lemma 5 in Ref.[24] shows this 
situation will arise if the normal vectors of the triangles vary 
by more than k/2 . A successful algorithm must guard 
against <t> inCircle becoming inconsistent. 

3.3. A New Surface Meshing Algorithm 

These outstanding ambiguities and computational expense 
motivated a search for a more straightforward algorithm. Our 
method makes two fundamental changes. 

• To avoid the ambiguity associated with the definition of 
^inCircle* we do not attempt a Delaunay triangulation of 
the curved surface. Instead we seek a tri angulation which is 
everywhere locally maxmin. While this may seem to consti- 
tute a dramatic relaxation, recall that the goal is a practical 
algorithm, and we have no reason to preler strictly Delaunay 
output triangulations. In addition, since one property of a 
Delaunay triangulation is that it is maxmin, this choice is 
worth investigating. 

• When ail angles of a triangle are acute, the circumcenter 
falls within the triangle. Ownership of the new site can then 
be uniquely assigned to this triangle. However, when a trian- 
gle is obtuse, ownership can become less clear. Therefore we 
make a simple choice: When a triangle is obtuse, we insert 
the new site at the centroid of the other triangle sharing the 
edge opposite the obtuse angle. Figure 4. illustrates this 
insertion rule. If tj is an obtuse triangle, then the new sites are 
added at the centroid of t opp . 

While admittedly ad-hoc , this site insertion strategy is not as 
arbitrary as it may initially seem. As a particular angle of tj 
opens up in the transition from acute to obtuse, the circum- 
center will pass from within t t to within t opp . The centroid of 
t may therefore be thought of as an approximate circum- 
center of tj. Denoting the radius of the circumcircle of tj as /?, 



tj is acute tj is acute 

Figure 4. Sites are added at the centroid of a triangle if the trian- 
gle is acute. Otherwise they are added at the centroid of the tri- 
angle opposite the obtuse angle. 
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it can be shown that if one chooses an approximate circum- 
center within a distance//? of the true circumcenter, then an 
angle bound of \ _ f 

a<asin( ^ ) (4) 

may still be achieved. Ref. [24] contends that even using 
approximate circumcenters, angle bounds of 30 may be 
achieved in practice. 

With these changes, the manifold triangulation algorithm 
becomes: 

Algorithm M: 

Quality Manifold triangulation of a CAD volume . 

Input: Underlying CAD volume, P, with initial triangulation 

S Cl and input vertices, V in . 

Target angle, a. 

Output: CX9iV out ) with all angles > a. 

Initialize: Compute constrained locally maxmin triangulation 
CX9^V iri )y using cad_edge entities as constraints. 
Build minimum angle priority queue, PQ mtn with t P Q 
denoting triangle at head ol queue, having min angle 

- ® p Q 

1. Apply <P Encroached to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CX9&V), Update PQ mnr 

1 

2. While (Q P q < a){ 

2. a If {tpQ is not obtuse) 

Assign t := t P g t with circumcenter p. 

Else Assign t:= t opp with centroid p. 

2.b If (p encroaches any constraining edge, e ) 

2.c Split constrained edge. Update CX9&V), Update PQ mnr 

2. d Else Split triangle /: 

Add p to V. Update CX^V), Update PQ mi(V 

} 

3. Output CXM VO. 

Note that when the algorithm terminates, it recovers the prop- 
erties of quality and size-optimality cited after the presenta- 
tion of Alg.Q in §3.1. The bound on output size, however, 
depends on the site insertion strategy always inserting new 
vertices within the circumcircle of 1 P q}^ and the modified 
strategy will not always guarantee this, thus the algorithm 
sacrifices strict proof of this property. 

The new algorithm requires initialization with a transforma- 
tion of the part’s original constrained manifold triangulation 
Sc to one which is everywhere locally maxmin. It we assume 
the existence of a triangulation predicate O xn which can be 
enforced for every pair of triangles sharing an edge in the tri- 
angulation (similar to the application of d > inCircle )> then 
this initialization may be performed with edge sweeps fol- 
lowed by edge-swapping when a violation is encountered. 
While the possibility of multiple sweeps makes this a seem- 
ingly inefficient approach, we recall that the initial complex- 
ity of S c is only proportional to the number of cad_edges in 
the geometric structure. Thus this simplistic approach is not a 
problem. 

Site insertion proceeds in a manner comparable to the PSLG 
algorithm. A convenient implementation makes use of the 
incremental insertion strategy of Ref.[25]. However, since the 
triangulation is no longer Delaunay, both forward and reverse 


propagation is necessary after edge swaps. Ref.[25] gives a 
modified point placement strategy for splitting encroached 
edges. Our implementation of Alg.M adopts this strategy 
without modification. 

3.4. Triangulation Predicates 

Algorithm M rests on two new triangulation predicates. One 
test a triangle for an obtuse angle, <t\ )/7 (0 and the other, 
<t> X/V (e) , tests if the edge, e , shared by any two triangles max- 
imizes the minimum angle in either triangle. This is accom- 
plished by comparison with a swapped edge, e\ which 
connects the opposite vertices of the two triangles. Our 
approach hinges on the hope of evaluating these predicates 
robustly in physical space. 

Both of these predicates can be formulated with direct mea- 
surement of the angles in a mesh. Since a triangle has three 
points, it uniquely define a plane in three dimensions. Angle 
measurement in a plane is unambiguous - despite the fact 
that the triangles form a discrete manifold which is lower 
dimensional than the surrounding space. Numerically, we 
recall that since vertex geometry returned by the CAD engine 
is considered exact, and so the error bound on the input is 
identically zero. 

<J>„, 7 (/) is applied by a logical or accumulation of 
i ■ j+ | , where j is a cyclic index running over the verti- 
ces of f, j € { 1 , 2, 3 } . This angle predicate is formed by com- 
parison of the square of the edge opposite j to that of a 
hypothetical edge formed by assuming that the edges of t 
incident on j form a right angle in the plane of t. 

Q> oh Vm) = v $oh( Zm ) v (5) 

where 

T if |*' y -v > . l | 2 + |v y .,-vJ i >|v / . l -v j _,p (6) 
F otherwise 

Notice that Eq.(6) requires only subtraction and multiplica- 
tion of data which is known exactly. Thus the computational 
requirements for exact evaluation are the same as for evalua- 
tion Of the dot product for <1 'encroached in §3.1. 


Evaluation of d>x/v(e ) requires a more direct method of angle 
measurement. Figure 5 shows a construction for this mea- 
surement. Recalling that sin( ) is monotone over the interval 





ij.y+i" 
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[0, tc/ 2] the construction in the figure shows that the differ- 
ence of the unit vectors of the edges incident upon any ver- 
tex, / is sufficient to define a vector d whose magnitude 
varies monotonically with the angle formed by the edges 
incident upon/ As in the preceding two predicates, computa- 
tionally, it is sufficient to evaluate only the square of this 
magnitude. 

The computational requirements for this predicate are some- 
what less simple. The computation ot unit vectors requires 
robust computation of the inverse of a vector magnitude 
\ /\V\ . Thus, exact evaluation ot <J y xN( e ) requires software 
arithmetic, like the packages in Refs. [26], [27] or [28]. 
Although this need must be viewed as a drawback, we note 
that it is confined to a single predicate, and to a relatively 
simple expression. Moreover, since the input geometry is 
exact, exact computation remains a feasible strategy. In the 
preliminary results shown in §4, all computation was per- 
formed using only double-precision floating- point hardware, 
and the option of software arithmetic was not pursued. 

3.5. Edge Refinement 

Algorithm M drives small angles out of the evolving triangu- 
lation. In addition we wish to satisfy both chord-height and 
length-scale criterion. After initially creating a triangulation 
free of small angles we apply an edge-refinement procedure 
to enforce these two requirements. Since the algorithm is the 
same for either of these criterion, Algorithm E considers a 
generic edge-based scalar y(e). In our implementation chord- 
height is defined as the square of the distance irom the mid- 
dle of an edge to the corresponding location on the actual sur- 
face of the model (provided through CAPRI by the geometry 
engine). 

Algorithm E: Edge Refinement of Manifold Triangulation. 

Input: Underlying CAD volume, P, with current triangulation, 

CXt 7\&V), and vertex set, V. 

Edge criteria y. 

Output: CX9iV out ) with all edges satisfying y(e) < y. 

Initialize: Build priority queue, PQy with c P q denoting edge at 

head of queue, having y{e P q). 

1. Apply ^Encroached to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CX%V), Update PQy 

) 

2. While (yie PQ ) >Y)( 

2. a Let p be midpoint ot e . 

2.b If (p encroaches any constraining edge, e) 

2.c Split constrained edge. Update CUiA/fVT Update PQy 

2. d Else Split edge e : 

Add p to e. Update CXA&VO, Update PQy 

} 

3. Output CXWY)- 

Assuming that y{e) is a static criterion, like chord-height tol- 
erance or length-scale, Alg.M can then be re-applied to 
CX'T&V) to remove any small angles created during edge 
refinement and swapping. 

3.6. Length-Scale Transport 

Quad/Octrees have been used extensively in unstructured 


mesh generation both for proximity testing and the establish- 
men! of local length scales. 129 ’ 30,31 ’ 32 ’” 1 The Cartesian mesh 
generation system of Ref.[2] provides just such a mechanism. 
Through the use of compact unstructured data structures, this 
particular Cartesian method preserves the ability to direction- 
ally refine Cartesian cells. 1341 Mesh generation schemes such 
as this produce geometry-refined meshes with smaller cells 
near smaller features ot the geometry. In addition, since it is 
an embedded-boundary approach, the Cartesian mesh gener- 
ator can outputs a set of cut-cells, C, which envelope the sur- 
face of the geometry. The triangulation output from Alg.M 
and chord-height refinement provides sufficient input for 
Cartesian mesh generation. 

Figure 6 demonstrates the geometry-adapted Cartesian 
mesh’s ability to transport length-scales between disjoint 
components. Since our intent is to produce a method which is 
fully automated, this property is essential. 

The set of N c cut-cells, C, is inserted into an Alternating Dig- 
ital Tree (ADT) [35] so that any edge, e. in CX9&V) can locate 
the subset of C which intersect it with a log(N c ) time bound. 
The ADT returns a short list of cells, which intersect e, 
C, = enC . Since these cells may have been created through 
anisotropic subdivision, we can compute an average size, h, 
for the intersected cells in each coordinate direction, 
je { 1, 2, 3} which defines a \ector length scale function for 
every edge in the manifold triangulation. 

(7) 

c ft = 1 


h{e) may then be used to compute the normalized excess 
length , lj ,of the edge in each coordinate direction. 



Figure 6. Isotropically refined Cartesian mesh around multi-ele- 
ment airfoil using method of Ref. [2]. Since the mesh is geo- 
metrically adapted, it communicates length scales between 
disjoint geometric entities in close proximity. 
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Figure 7. isotropic (top) and anisotropic (bottom) sets ot cut- 
cells from Cartesian mesh generator used to specify length 
scales for surface triangulation of manufacturing die example. 



e: -hXe) 


( 8 ) 


Z* = Max{l\, / 2 , h) (9) 

The scalar /* is computed for each edge and inserted into the 
priority queue PQy in Alg.E for length-scale based edge 
refinement. 


Figure 7 displays the set of cut-cells from isotropic and 
anisotropically refined Cartesian meshes for the manufactur- 
ing die example from §2.2 (Fig. 2). This set of cells defines 
the vector length scale function for the surface via Eq.(7). 


3.7. An Alternative Insertion Strategy 


The cut-cells in Fig. 7 are interesting in that they suggest an 
alternative to the incremental mesh improvement algorithm 
of the preceding sections. Rather than incremental site inser- 
tions, one can imagine an approach which simply includes 



Figure 8. An alternative site insertion strategy based upon pierce 
locations of the “most normal” edges from cut Cartesian cells 
and cad_edge/Cartesian cell face pierce points. 


the intersection points of the “most normal” Cartesian edges 
with the surface triangulation. In order to preserve the 
cad_edge features of the geometry, we augment this set with 
the pierce locations of cad_edges with the faces of the cut- 
cells. Figure 8 demonstrates the construction of a vertex set 
in this manner. 

We explore this alternative using the flap geometry and cut- 
cell set for the flap geometry shown in Figure 9. The surface 
triangulation formed by enforcing i>xN on this point set is 
shown in Figure 10. In examination ot this figure several 
things are apparent. Despite the fact that the flap is swept, 
twisted and inclined, the anisotropic Cartesian cell subdivi- 
sion correctly introduces high-aspect ratio cells along the 
flap, with vertices clustered at the leading and trailing edges 
as a result of the geometry adapted Cartesian mesh. Never- 
theless, the triangulation is somewhat irregular, both as a 
result of swaps enforced by the <&xn predicate and at loca- 
tions corresponding to the interfaces between adaptation lev- 
els of the cut Cartesian cells. Other triangulation predicates 
(Minmax, for example) may alleviate some of this irregular- 
ity,^ 7 ! 

A truly successful implementation of this insertion strategy 
would require explicit mesh smoothing passes to mitigate 
these irregularities. In addition, in examining the triangula- 




Figure 10. Surface triangulation formed using sites of “most 
normal” Cartesian edges from Fig.9 and cad_edge pierce 
locations from Fig. 8. 
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tion at the leading edge (close-up, Fig. 10) we can discern a 
region of irregularity associated with the switch over ol 
pierce locations from horizontal Cartesian edge-pierces to 
vertical edge-pierces. Mesh smoothing would be beneficial in 
this area as well. 

One very attractive property of this approach is that it com- 
pletely decouples the final triangulation from the description 
of the geometry. The final vertex set is not a function of the 
initial vertex set. Therefore, any “history is removed from 
the triangulation. The surface triangulation in Fig. 10 should 
be considered preliminary, and although it has some attrac- 
tive properties this approach was not aggressively pursued. 

4. Results and Discussion 

This section presents example meshes on several CAD parts 
from a variety of sources. All the example parts were read in 
their native CAD tile format using the CAPRI API without 
special treatment. The investigations focus on examination of 
issues raised in the presentation of the triangulation algo- 
rithm in §3.3, the edge refinement strategy from §3.5 and the 
length scale transport from §3.6. 

4.1, Minimum Angle Bound 

In §3.3, Algorithm M was presented without firm proof of 
termination. Moreover, the discussion noted that the modified 
site insertion strategy for obtuse triangles violates one of the 
assumptions that establishes a bound on the output size of the 
mesh in the PSLG method. It is therefore necessary to dem- 
onstrate the performance of the Alg.M to show that it both 
terminates and produces meshes with an economy similar to 
that of the PSLG method upon which it is based. 

Figure 1 1 contains a histogram of the evolution ol the small- 
est angle in the mesh as Alg.M proceeds on the manufacture 


30 



5 


0 i .. 1 — — — 

0 2000 4000 6000 8000 1 0000 

Number of Triangles in Manifold 

Figure 11. Histogram of minimum angle during mesh evolution 
using Alg.M (without the edge refinement of §3.5) for the manu- 
facturing die example presented earlier. In this example, a mesh 
with a minimum angle of 25° would contain 2606 triangles. 



Figure 12. Quality manifold triangulation for manufacturing die 
example generated with Alg.M. in §3.3. Mesh improvement 
terminated after generating 2606 triangles when the mini- 
mum angle in the triangulation reached 25°. Chord-height 
and length-scale refinement not used. 

ing die example problem used in earlier illustrations. While 
this curve is far from monotone, it clearly displays the steady 
improvement of the minimum angle in the mesh as the site 
algorithm proceeds. The steep initial rise indicates rapid 
annihilation of extremely small angles in the mesh, and the 
mesh achieves a minimum angle of almost 29° by the end of 
the histogram. 

The dashed line at 25° highlights the first time that all angles 
in the mesh exceeded this value. Tracing this value on the 
abscissa shows that setting the angle bound, a, to 25 c will 
cause Alg.M to terminate after generating 2606 triangles. 
Figure 12 shows the resulting triangulation. As discussed in 
§3.1 and §3.3, the presence of an angle bound ensures that 
small features are surrounded by proportionally small trian- 
gles (see the inset frames in Fig. 12), and that the mesh length 
scale varies smoothly over the part. The smallest angle in the 
mesh is 25.02°, which corresponds to a maximum aspect 
ratio of between 2.36 and 4.7 (by Eq.(3)). 

While the histogram in Figure 1 1 shows a steady increase in 
the minimum angle, there is an irregular array of downward 
spikes in the profile. In the presentation of the PSLG algo- 
rithm, Ref.23 included a similar histogram, and noted the 
same characteristic. Consider the two triangles t\ and 
shown at the left in Figure 13, Assume that 0j is the smallest 



Figure 13. Mechanism responsible for downward spikes in histo- 
gram of minimum angle as Alg.M proceeds. If 0j is initially the 
smallest angle in the mesh, angles 0 3 and 0 4 may be smaller 
after insertion of p and enforcement of Q>xn by edge-swapping. 


9 



- AlAA 99-0776 - 


angle in the mesh lying in triangle Furthermore, assume 
that ffs face neighbor, t 2 , has a small angle 0 2 opposite the 
shared edge which is only slightly larger than Bp Site p will 
get added at tf s circumcenter improving 0! to 2 0j. After 
application of the maxmin predicate <t>xyv on l ^e shared edge, 
the swapped configuration at the right of Fig. 13 may occur. 
This configuration includes two new triangles with a mini- 
mum angles 03 and 04 either of which may now be the small- 
est angle in the mesh and may actually be smaller than the 
original angle 0j. 

With this behavior understood, Figures 14 and 15 present 
angle histograms and example meshes for a more compli- 
cated set of parts. Alg.M was run on CAD parts for the mam 
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Figure 14. Angle histograms for Alg.M on the main element of a 
transport wing(upper), and on the main flap (lower). An angle 
bound of 25° would produce tri angulations with 20846 and 
15334 triangles on the main element and flap respectively. 



Figure 15. Bounded angle triangulations of main element of a 
transport aircraft wing and flap generated by Alg.M in §3.3. 
Minimum angle 25°, 20846 triangles on wing, 15334 triangles 
on flap. 

element of a transport wing, and a flap element for the same 
wing. The main element consisted of 224 rational B-spline 
curves and 36 trimmed NURBS surfaces. The flap contained 
31 rational B-spline curves and 10 trimmed NURBS sur- 
faces. 

The crosshairs on the curves in Fig. 14 show that with a 25 
angle bound, Alg. M will produce 20846 triangles on the 
main element and 15334 triangles on the flap. Figure 15 dis- 
plays these triangulations. 

The histograms in Figure 14 bear close resemblance to the 
one presented for the simple die example shown in Fig. 1 1. 
All of these profiles are characterized by a sharp initial angle 
improvement and then a rolling-off as the minimum angle 
climbs above about 20°. All three profiles exceed 27°, but 
reaching 30° seems unlikely. 

The abscissa values in Fig. 14 indicate relatively large trian- 
gulations as compared with the die example. The size-opti- 
mality property cited in §3.1 and §3.3 implies that triangle 
size must vary smoothly between different sized constraints. 
In both of these examples the smallest constraint in the capri 
triangulation was a factor of 10^ smaller than the largest con- 
straining edge. In addition, practical experience with the 
algorithm indicates that larger initial triangulations, S 0 
(Eq.(2)) have a proportionally slower initial rise in their angle 
profiles. This seems reasonable since if the CAPRI triangula- 
tion is complex there may initially be many bad angles which 
need improvement. Often CAD parts are unnecessarily com- 
plex for reasons dating to the specific events in their creation. 
Our experience with CAD repair software indicates that 
CAPRI generally produces less complex (sometimes by an 
order of magnitude) initial triangulations if the parts have 
been processed by CAD repair tools. Since tolerance to poor 
CAD parts is one type of robustness that we seek in this 
research, neither part shown in Fig. 15 underwent such repair 
prior to creation of these triangulations. 
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Figure 16. Interior corner of manufacturing die example trian- 
gulation after imposition of a chord-height tolerance ot 
l.xlO' 4 normalized by the maximum outer dimension of the 
part. Edge refinement was performed using Alg.E. in §3.5 

With the minimum mesh angle restricted to 25°, triangles are 
again restricted to aspect ratio’s less than 4.7. Such an isotro- 
pic mesh is very inefficient at meshing features with curva- 
ture in only one dimension. Therefore, if one intends to 
produce meshes for viscous computation, a substantially 
smaller angle bound may be appropriate. 

4.2. Chord-Height Refinement 

Enforcement of the angle bound does not guarantee that the 
edges are refined when they are sufficiently far from the 
underlying surface. Alg.E in §3.5 presented an edge refine- 
ment strategy for automatically breaking edges whose mid* 
points which are far from the geometry. Figure 16 shows an 
interior corner of the die example after the imposition of a 
chord-height tolerance of 10" 4 times L, where L is the maxi- 
mum outer dimension of the part. Away from the curved inte- 
rior comers, the triangulation remains unchanged since the 
faces of the die are planar. 

4.3. Length Scale Refinement 

Once the initial triangulations of Figs. 12 and 15 have been 
supplemented with the addition of sites to match the chord- 
height bound, the triangulations accurately describe the 
geometry of the part to be meshed. These triangulations then 
need refinement to reflect the length-scales brought in by the 
cut Cartesian cells. 

The set of anisotropically refined cut-cells from Fig. 7 was 
used to define a vector length scale function for each edge in 
the triangulation using eq.(7). Edges were processed with 
Alg.E using a priority queue based on the maximum compo- 
nent of the normalized excess length function from eqs.(8) 
and (9), PQ { *. Figure 17 displays the resulting length-scale 
adapted mesh for the die example. This manifold triangula- 
tion shown contains 17384 triangles and processing of the 
PQ { * was halted after achieving a value of /*= 1.0. In addi- 
tion to the /*, based edge refinement, the minimum angle 
bound was set to 15° and a chord-height tolerance of 
U10“ 4 L. In comparing this triangulation with its earlier 
forms (in Figs. 12 and 16) one sees substantial refinement ot 
the triangulation in response to the length-scale information. 
However, such a comparison also shows that application of 
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Figure 17. Top, bottom and 3- view of manufacturing die 
example geometry after length-scale refinement using 
the anisotropic Cartesian cut-cells from Fig.7. 17384 tri- 
angles, a = 15°, maximum chord-height set to 1 .xlO' 4 L. 

Alg.E has apparently failed to convey the anisotropy of the 
cut-cells to the final triangulation (c.f Fig.7). The reasons 
behind this shortcoming are obvious. 

(1) Our edge-refinement uses the very simplistic strategy 
of simply breaking edges at their midpoint, without con- 
sulting Ji(e) for information regarding desired anisot- 
ropy prior to point placement. 

(2) The triangulation predicate 6*^ lights anisotropy after 
every insertion, since it forces edge swaps in order to 
maximize the minimum angle. 

In order to use this approach for creating acceptable stretched 
triangulations, then the triangulation predicate and point 
placement strategies need to be modified. One interesting 
approach would be to modify i>xN with h(e) by applying a 
scaling during evaluation. In this way, the edge-swapping 
will be cognizant of the local stretching information. Such an 
approach is conceptually quite similar to the locally-mapped 
Delaunay methods sometimes used in 2-D Ref.[39] 
examined the performance of alternative triangulation predi- 
cates combined with modified site insertion rules and demon- 
strated that they can be effective in generating smooth 
stretched triangulations. Both of these approaches merit 
investigation since they fit well within the framework of the 
present algorithm, especially since the stretching information 
in h(e) is currently not utilized. 

Figure 18 contains length-scale adapted triangulations for the 
main wing element and flap geometries shown earlier in 
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Figure 18. Transport aircraft wing and flap geometries 
from Fig. 15 after length-scale matching as described in 
§3.6. The wing and flap have 34720 and 26640 triangles 
respectively. Inset frame shows matching length-scales at 
the flap gap. Minimum angle 15.03°. Chord-height toler- 
ance set to 1 x 1 0"^ of the span. 

Fig. 15. These triangulations were again created using a mini- 
mum angle bound set to 15 which resulted in a smallest 
angle of 15.03°. A chord-height tolerance of 1 a 1 0~ L , was 
used, however, this value is somewhat lax for a serious com- 
putation. Length-scale refinement was stopped at /*= L0. 
The triangulation of the main element has 34720 triangles 
while the flap has 26640. In comparing these triangulations 
to those in Fig. 15 notice that the flap leading-edge length 
scales have been transported through the flap gap to the flat 
surfaces of the main element. The small lrame inset in Fig. 
18 is included to show this detail. 

5. Conclusions and Future Work 
5.1. Conclusions 

This work focused on the direct use of CAD geometry for the 
automatic generation of closed manifold surface triangula- 
tions. Our approach used the CAPRI API to provide a mod- 
eler independent method of interacting directly with the CAD 
system’s geometry engine. This approach avoids data transla- 
tion which can deplete a model of topological information, 
and avoids the consistency conflicts which can occur when 
different geometry engines attempt to infer topology from a 
parts geometry. In addition, it permits the generation of new 
vertices using the same constructors which were used to cre- 
ate nearby geometry already in the part, thus avoiding build- 
ing an inconsistent model. 

A second novel aspect of this work was the development of a 
new curved surface meshing technique based upon the use of 
approximate circumcenters for site insertion. The method 
maintains a locally maxmin triangulation and produces a 
bounded aspect ratio manifold triangulation. The algorithm 
was demonstrated on CAD parts of varying complexity and 
reliably produced triangulations with minimum angles in 
excess of 27°. 


An edge refinement algorithm was also presented which was 
used to ensure that edges in the triangulation match the 
underlying part model to within a specified tolerance. Edge 
refinement was also used in conjunction with an anisotropi- 
cally refined Cartesian mesh to ensure that triangles on dis- 
joint surfaces in close proximity maintain matching length- 
scales. 

An alternative site insertion strategy was examined in which 
the final vertex set was taken from the “most normal pierce 
points of a geometry-adapted Cartesian mesh with an early 
triangulation. This approach remains attractive since it fully 
decouples the final triangulation from any intermediate mod- 
els used to initially describe the surface. 

5.2. Future Research 

• Despite the fact that the cut Cartesian cells define a vector 
length scale function in Eq.(8), our simple (edge subdivision) 
insertion strategy does not make use of this information. 
Since the adapted Cartesian mesh was subdivided in response 
to surface curvature, hj(e) should be able to provide more 
guidance concerning the placement of new sites. Since we 
intend to produce surface triangulations that include high 
aspect ratio triangles finer control of site placement must is a 
necessity. Two strategies for including this information were 
cited in §4.3, as well as alternative triangulation predicates. 

• Although Alg.M. produces meshes with generally smooth 
length-scale variation. There is occasionally a discernible 
irregularity in the triangulations. This behavior becomes 
more pronounced when higher angle bounds are specified. 
Similar behavior has been noted in the PSLG algorithm [23], 
although in 2-D the behavior seems less pronounced. One 
possible source for this is the abrupt change in insertion loca- 
tion from circumcenter to centroid (of t 0f? p) if an obtuse trian- 
gle is encountered. Alternative strategies should be 
investigated since this behavior can degrade the efficiency of 
the triangulations. 

• The initial triangulation returned by CAPRI depends 
strongly upon the “history” of the CAD part. Such effects are 
evident in the triangulation of the main wing element in 
Fig. 15, for example, where the outboard portion of the flap 
cut-out is excessively resolved due to relics of the part’s cre- 
ation process. The only effective strategy we’ve found for 
avoiding this is to “clean” the geometry using CAD repair 
software. Alternative approaches should be investigated. 
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